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/  ABSTRACT 

/ 

'  /  Splitting  methods  provide  efficient  tools  for  solving  linear  and 

nonlinear  time  dependent  problems  modelled  by  partial  differential 
equations.  In  this  report  we  discuss  the  numerical  solution  of  the  Navier- 
Stokes  equations  for  incompressible  viscous  fluids  by  such  methods.  The 
splitting  permits  decoupling  the  two  main  difficulties  in  the  problem,  namely 
the  nonlinearity  and  the  incompressibility.  Actually  these  splitting  methods 
have  a  broad  range  of  applicability  and  can  be  applied  for  example,  to  the 
solution  of  eigenvalue  problems.  ’  ’! 
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SIGNIFICANCE  AND  EXPLANATION 


The  numerical  solution  of  the  Navier-Stokes  equations  for  Incompressible 
viscous  fluids  is  a  complicated  problem,  particularly  for  sufficiently  large 
Reynold  numbers  where  strong  nonlinear  effects  are  present.  Usii  g  an 
appropriate  operator  splitting  method,  we  show  in  this  report  that  it  is 
possible  to  decouple  the  two  main  difficulties  in  the  problem,  namely  the 
nonlinearlity  and  the  incompressibility  condition.  This  reduces  the  solution 
of  the  time  dependent  problem  to  a  sequence  of  simpler  stationary  problems, 
which  can  be  solved  in  many  cases  by  standard  procedures  such  as  conjugate 
gradient  algorithms,  simple  finite  element  approximations,  etc...  . 

It  is  also  shown  in  this  report  that  these  splitting  methods  can  be 
applied  to  eigenvalue  calculations  and  in  fact,  they  have  been  applied 
elsewhere  to  the  solution  of  nonlinear  eigenvalue  problems  such  as  the  Hartree 
equation  in  Quantum  Physics. 

Numerical  results  show  the  possibility  of  the  solution  methods  discussed 
in  this  report  when  it  comes  to  the  simulation  of  flows  closely  related  to 
practical  problems. 

On  the  basis  of  these  numerical  results  and  accompanying  experiments  it 
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SPLITTING  METHODS  FOR  THE  NUMERICAL  SOLUTION  OF  THE 
INCOMPRESSIBLE  NAVIER-STOKES  EQUATIONS 


R.  Glowinski 


1.  Introduction.  Synopsis 


Solving  by  numerical  methods  the  Navier-Stokes  equations,  in 
order  to  simulate  unsteady  flows  of  either  incompressible  or  compres¬ 
sible  viscous  fluids ,  is  still  a  challenging  problem.  This  important 
problem  has  motivated  the  work  of  many  scientists  (see,  e.g.  TEMAM  [1  ]  , 
GIRAULT-RAVIART  [2]  ,  RAUTMANN  [3]  ,  THOMASSET  [4]  ,  GLOWINSKI  [5 3  for 
references).  Concentrating  on  the  incompressible  case  we  would  like  to 
show  in  this  paper  that  operator  splitting  methods,  like  those  advocated 
by  Professor  G.I.  Marchuk  in  [6 3  ,  provide  quite  efficient  numerical  sche¬ 
mes  for  solving  the  time  dependent  Navier-Stokes  equations.  The  content 
of  the  paper  is  as  follows  ; 


In  Section  2,  we  describe  and  comment  the  Navier-Stokes  equations 
modelling  unsteady  flows  for  incompressible  viscous  fluids.  In  Section  3 
we  discuss  some  general  schemes  using  operator  splitting  and  apply 
them  to  the  Navier-Stokes  equations  in  order  to  decouple  incompres¬ 
sibility  and  nonlinearity  (in  fact  we  could  not  resist  concluding 
that  section  by  showing  that  the  same  principles  also  apply  to 
eigenvalue  calculations).  In  Sections  4  and  5,  we  discuss  the  speci¬ 
fic  treatment  of  the  nonlinearity  and  of  the  incompressibility,  res¬ 
pectively,  Finite  Element  Approximations  are  discussed  in  Section  6, 
and  finally  we  show  in  Section  7  the  results  of  several  numerical  ex¬ 
periments  designed  to  test  the  methods  previously  discussed  in  the 
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2.  Formulation  of  the  unsteady  Navier-Stokes  equations  for  incompres¬ 


sible  viscous  fluids. 


Let  us  consider  a  newtonian  viscous  and  incompressible  fluid. 

If  ft  and  T  denote  the  region  of  the  flow  and  its  boundary,  respec¬ 
tively,  then  this  flow  is  governed  by  the  Navier-Stokes  equations 


(2.1)  37  "  vAu  +  (u.V)u  +  Vp  -  f  in  (1, 


(2.2) 


V.u  -0  in  ft  (incompressibility  condition ). 


In  (2.1),  (2.2)  : 


N 

(i)  u  ■  {u. is  the  flow  velocity, 
1  1“  1 


(ii)  p  is  the  pressure. 


(iii)  v  is  the  viscosity  of  the  fluid  (in  normalized  units  we  have 
v  -  1/Re,  where  Re  is  the  Reynold's  number), 


(iv)  f  is  a  density  of  external  forces. 


In  (2.1),  (u.V)u  is  a  symbolic  notation  for  the  nonlinear  (vector) 


term  : 


N  3u.  N 

{ ^  u;  3”}.  i  • 

j-i  J  j  1-1 


Boundary  conditions  have  to  be  added  ;  for  example,  in  the  case  of 
the  airfoil  A  of  Fig.  2.1,  below,  we  have  (since  the  fluid  is  viscous) 
the  following  adherence  condition 


(2.3) 


u  ■  0  on  3A  ■  T.  . 

^  A 


TOWS1 


.V.’iU 


I? 


SPLITTING  METHODS  FOR  THE  NUMERICAL  SOLUTION  OF  THE 
INCOMPRESSIBLE  NAVI ER-STOKES  EQUATIONS 


R.  Glowinski 


1 .  Introduction.  Synopsis 

Solving  by  numerical  methods  the  Navier-Stokes  equations,  in 
order  to  simulate  unsteady  flows  of  either  incompressible  or  compres¬ 
sible  viscous  fluids,  is  still  a  challenging  problem.  This  important 
problem  has  motivated  the  work  of  many  scientists  (see,  e.g.  TEMAM  [13, 
GIRAULT-RAVIART  [2  3  ,  RAUTMANN  [3  3,  THOMAS  SET  [4  3,  GLOWINSKI  [5  3  for 
references).  Concentrating  on  the  incompressible  case  we  would  like  to 
show  in  this  paper  that  operator  splitting  methods,  like  those  advocated 
by  Professor  G.I.  Marchuk  in  [63,  provide  quite  efficient  numerical  sche¬ 
mes  for  solving  the  time  dependent  Navier-Stokes  equations.  The  content 
of  the  paper  is  as  follows  : 
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In  Section  2,  we  describe  and  consent  the  Navier-Stokes  equations 
modelling  unsteady  flows  for  incompressible  viscous  fluids.  In  Section  3 
we  discuss  some  general  schemes  using  operator  splitting  and  apply 
them  to  the  Navier-Stokes  equations  in  order  to  decouple  incompres¬ 
sibility  and  nonlinearity  (in  fact  we  could  not  resist  concluding 
that  section  by  showing  that  the  same  principles  also  apply  to 
eigenvalue  calculations).  In  Sections  4  and  5,  we  discuss  the  speci¬ 
fic  treatment  of  the  nonlinearity  and  of  the  incompressibility,  res¬ 
pectively,  Finite  Element  Approximations  are  discussed  in  Section  6, 
and  finally  we  show  in  Section  7  the  results  of  several  numerical  ex¬ 
periments  designed  to  test  the  methods  previously  discussed  in  the 
paper . 
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2.  Formulation  of  the  unsteady  Navier-Stokes  equations  for  mcompres 


sible  viscous  fluids. 

Let  us  consider  a  newtonian  viscous  and  incompressible  fluid. 

If  ft  and  T  denote  the  region  of  the  flow  and  its  boundary,  respec* 
tively,  then  this  flow  is  governed  by  the  Navier-Stokes  equations 


(2.1)  -r~  -  vAu  +  (u.V)u  +  Vp  ■  f  in  ft, 

d  tl  ~  — 

(2.2)  V.u  -0  in  ft  (incompressibility  condition ). 

In  (2.1),  (2.2)  : 

N 

(i)  u  *  {u^}^  velocity, 

(ii)  p  is  the  pressure, 

(iii)  v  is  the  viscosity  of  the  fluid  (in  normalized  units  we  have 
v  -  1/Re,  where  Re  is  the  Reynold's  number), 

(iv)  f  is  a  density  of  external  forces. 

In  (2.1),  (u.V)u  is  a  symbolic  notation  for  the  nonlinear  (vector) 
term  : 

N  du.  N 

{I  u.  ■). 

j-1  3  i-' 

Boundary  conditions  have  to  be  added  ;  for  example,  in  the  case  of 
the  airfoil  A  of  Fig.  2.1,  below,  we  have  (since  the  fluid  is  viscous ) 
the  following  adherence  condition 


(2.3) 


u  -  0  on  9A  ■  T  . 
-  A 
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Typical  conditions  at  infinity  are 

(2.4)  u  •  uM  , 

where  ub  is  a  constant  vector  (with  regards  to  the  space  variables 
at  least) . 

Finally,  for  the  time  dependent  problem  (2.1),  (2.2),  an  initial  condition 
such  as 

(2.5)  u(x,0)  -  uq(x)  on  (2  , 

where  uq  is  given  (with  V.uo  *0),  is  usually  prescribed. 

Uith  the  methods  described  in  this  paper,  we  can  also  treat  boundary 
conditions  such  as 

3u 

(2.6)  u  -  go  on  rQ  ,  v^-np-gj  on  T ,  , 


for  the  numerical  implementation  of  (2.6)). 

In  two  dimensions,  it  may  be  convenient  to  formulate  the  Navier-Stokes 
equations  using  a  stream-function-vorticity  formulation  (see,  e.g., 
FORTIN-THOMASSET  [7]  ,  GIRAULT-RAVIART  [2  ]  ,  GLOWINSKI-PIRONNEAU  [8  ]  , 
REINHART  [9]  ,  KELLER-SCHREIBER  [10]  ,  GLOWINSKI-KELLER-RE INHART  [11]  , 
ROACHE  [12  ]  ) . 

Solving  the  above  Navier-Stokes  equations,  even  at  moderate 

3 

Reynold  number  (say  Re  -  10  ),  is  a  non  trivial  task  because  the 
following  difficulties 

(i)  It  is  a  nonlinear  system  of  partial  differential  equations, 

(ii)  We  have  the  inconpressibility  condition  (2.2). 

In  the  following  Section  3,  we  shall  see  how  time  discretization  by  opera¬ 
tor  splitting  methods  can  decouple  the  nonlinear  and  incompressibility 
difficulties. 

To  conclude  this  section,  let  us  mention  that  a  mathematical  ana¬ 
lysis  of  the  Navier-Stokes  equations  for  incompressible  viscous  fluids 
can  be  found  in  e.g.  LIONS  [13]  ,  LADYSHENSKAYA  [14]  ,  TEMAM  [1  ]  and 
TARTAR  [15]  . 


(3.2) 


A  ■  Aj  +  A2 


( non  trivial  meaning  that  A ^  and  A2  are  individually  simpler  than  A) . 
With  At  (>  0)  a  time  discretization  step,  let  us  define  several  sche¬ 
mes  taking  advantage  of  the  decomposition  (3.2)  : 


A.  A  Peaceman-Rachford  time  discretization  scheme. 
The  scheme  is  defined  as  follows 


(3.3) 


u  -  uo  , 


then  for  n  z  0,  with  un  known ,  we  compute  successively  un+1^2  and 
then  un+1  by 

,,  ,s  Un+1/2-  un  .  .  ,  n+l/2s  .  .  ,  nx  ,n+l/2 

(3.4)  - -  +  Aj(u  )  +  A2(u  )  -  f  , 


(3.5) 


un+1  -  un+1/2  .  n+1/2.  /  n+ 1 .  fn+l 

- ^72 - +  Al(u  )  +  A2(u  )  f 


In  (3.4),  (3.5)  u  denotes  an  approximation  of  u((n+ci)At) ,  and 
fn+a  -  f((n+a)At). 


B.  A  Douglas-Rachford  time  discretization  scheme 

It  is  the  variant  of  the  above  scheme  described  by 


(3.6) 


u  -  uo, 


t hen,  for  n  i  0,  define  fln+*  and  un+  *  from  un  by 


(3.7) 


~n+ In.  . 

u  -  u  .  .  ,«n+K  .  .  ,  nN  „n+l 
- ^ -  +  Aj(u  )  +  A2(u  )  -  f  , 


n+ 1  n  .  ,  , 

„  o^  u  -  u  .  .  ,An+ 1 s  .  .  ,  n+lx  „n+l 

(3-8)  - ^ -  +  Aj(Q  )  +  A2(u  )  -  f 


C.  A  three  stages  operator  splitting  scheme. 


Let  6  belongs  to  the  open  interval  (0,1/2)  ;  the  idea  behind  the 
scheme  is  to  split  the  time  interval  [nAt, (n+l)At  ]  in  three  subinter¬ 
vals,  as  shown  on  Figure  3.1, 


(n+8)At 


(n+1-8) At 


nAt 


¥r 


* - 1 - =► 

(n+l)At  1 


Figure  3. 1 . 


and  integrate  in  time  using  an  implicit  scheme  for  Aj  (resp.  an  explicit 
scheme  for  A2)  on  [nAt, (n+0)At 3  ,  then  switch  the  role  of  Aj  and  A2  on 
[(n+8)At,  (n+l-8)At  ]  ,  and  on  [(n+l-8)At,  (n+l)At  ]  do  like  on 
[nAt,  (n+8)At]  .  Using  these  principles,  we  obtain  the  following  scheme, 
some  forms  of  which  have  been  advocated  by  STRANG  [16]  ,  BEALE -MAJDA  [17] 
LEVEQUE  [18]  ,  LEVEQUE-OLIGER  [19]  (for  8  -  1/4)  s 


(3.9) 


o 

u  -  u0. 


then ,  for  n  k  0,  we  obtain  un+®  ,  un+1  ®,  un+1,  from  un,  as  follows 


(3.10) 


un+®-  u11  .  .  /..n+8x  .  »  /  nx  ,n+8 

-~eK — +  A1(U  >  +  A2  u  >  "  f 


n+1-8  n+8 
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n+l  n+1-8 


n+1-8.  ,n+l-8 


(3.12) 


u  -  u 


8At 


*  /  a+l\  ^  *  /  n+1-8.  rn+l 
+  Aj(u  )  +  A2(u  )  -  f 


The  convergence  of  (3.3)-(3.5)  and  (3.6)-(3.8)  has  been  proved  in 
LIONS-MERCIER  [20]  (see  also  GODLEWSKY  [21  ]  )  under  quite  general 
monotonicity  assumptions  on  Aj  and  A2  ;  it  is  very  likely  that  the 
methods  used  in  [20  ]  to  prove  the  convergence  of  (3.3)-(3.5)  and  (3.6)- 
(3.8)  still  apply  to  (3 . 9) — (3 . 12) . 

3.2.  Convergence  and  stability  properties  of  the  basic  schemes 
Following  the  approach  used  in  [6  ]  ,  we  shall  consider  for  sim¬ 
plicity  the  case  where  H  ■  f  ■  0,  u  e  1R1*,  A  is  an  N  x  N,  symmetric 

and  positive  definite  matrix  and  where 


(3.13)  Aj  ■  aA,  A2  -  BA  ,  with  cl  +  B  ■  1,0  <  a,B  <  1  . 
In  that  case  the  solution  of  (3.1)  is  clearly  given  by 


We  have,  from  (3.4),  (3.5),  (3.13), 


(3.15)  un+1-  (1+  6  A)"1  (I-  ot  A)  (1+ 


~  At  /r  o  At  A.  n 
Ot  -y  A)  (1-6  -y  A)u 


Using  a  vector  basis  consisting  of  eigenvectors  of  A,  we  have  from 
(3.15),  with  obvious  notation, 


(3.16) 


n+1 

u. 

l 


(1-  a  ~  A.)(l-  6  y-  A.) 

(1+  a  ^  A.)(l+  6  ^  Ui 


where  A.  (>  0, 

i 

that  Aj  s  A2  i, 
by 


Vi  -  1,...,N)  is  the  ith  eigenvalue  of  A  ;  we  suppose 
. ..  An  .  Consider  now  the  rational  function  Rj  defined 


(3.17) 


Rj  (x) 


a-  \  xxi- 1  x) 

(1+  |  x)(l+  |  x) 


we  observe  that  |Rj(x)|  <1  Vx  >  0,  implying,  in  that  simple  case,  the 
unconditional  stability  of  scheme  (3.3)-(3.5).  Since 


(3.18)  lim  R. (x)  -  1, 

X  -H-  oo  1 

we  observe  that  for  stiff  problems,  i.e.  problems  such  that  A^/Aj  »  1, 
scheme  (3.3)-(3.5)  is  not  very  good  to  damp,  simultaneously,  the  compo¬ 
nents  of  u11  associated  to  the  large  and  to  the  small  eigenvalues  of  a  ; 
from  this  observation,  we  can  expect  that  scheme  (3.3)-(3.5)  is  not 
well  suited  to  capture  the  steady  Btate  solutions  of  stiff  problems 
(like  those  obtained  from  the  discretization  of  partial  differential 
equations)  ;  this  has  been  confirmed  by  numerical  experiments. 

Since 

2 

(3.19)  e"x  -  1  -  x  +  +  x2e(x) 

and,  from  (3.17), 


f-.'- 


K*S 
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(3.20) 


x  2 

Rj  (x)  -  i-x  +  y~  +  x  n(*)  * 


with  lim  C(x) 
x  -*■  0 


lim  ri(x)  ■  0,  we  have  that  scheme  (4.3)-(4.5)  is 


aeoond  order  accurate  in  the  simple  case  that  we  have  considered. 

We  observe,  from  (3.15)  that  if  one  takes  a  ■  8  ■  1/2,  then  the  two 
linear  systems  which  have  to  be  solved,  at  each  full  step,  are  in  fact 
associated  to  the  same  matrix  1+4^-  A. 

Analysis  of  scheme  (3.6)-(3.8)  :  We  have  this  time 

(3.21)  un+1  -  (1+  aAtA)-1(I+6AtA)"1(I+  aB|At|2  A2)un  , 


which  implies 


j  1+  a6 | At | 2  X2 


(3.22)  u. 


n+ 


u? 


1  (1+  aAt  Xi)(l+  BAt  JL)  1 

If  we  define  now  R 2  by 

(3.23)  R,(x)  -  1  *  -B-*2 -  , 

L  (1+ax) (1+Bx) 

we  observe  that  0  <  R^(x)  <  1  ¥x  >  0,  implying  in  turn,  in  the  case 
under  discussion,  the  unconditional  stability  of  scheme  (3.6)-(3.8). 
Since  we  have  again 


(3.24)  lim  R,(x)  -  1, 
x  +» 

scheme  (3.6)-(3.8)  may  behave  poorly  for  stiff  problems  and  be  not  too 
efficient  for  capturing  steady  state  solutions.  About  the  accuracy  of 
scheme  (3.6)-(3.8),  we  should  essily  prove  that 

(3.25)  RjOO  “  1-x  +  x2+  x2  rj(x), 

with  lim  n(x)  ■  0,  implying  if  we  compare  to  (3.19),  that  scheme 
x  -*■  0 

(3.6)-(3.8)  is  only  first  order  accurate. 


*.\»  ,•*,*  /  V*  •'  "  ■  *  .  •  •  -  *  *.  *.  .  %  A  N,  *,  •*  %  .*»  A  A  A  „%*».,%•  A 


The  good  choice  £or  a  end  0  is  egein  a  ■  6  ■  1/2. 


Analysis  of  scheme  (3.9)-(3. 12) t  We  have  (with  6'  ■  1-20) 


(3.26)  « 


n+1-(I+  CX0At  A)_1(I-00At  A)(I+00'At  A)"1  (I-a0'AtA) 
x(l+a6  At  A)"1 (1-00  At  A)un  , 


which  implies 

(1-00 At  X.)2(l-a0'At  X.)  w 

/o  o  n+i  x  in 

(3.27)  u*  m  x  u*  . 

1  (l+a0At  Xi)z(l+00'At  Xt)  1 
Consider  now  the  rational  function  R^  defined  by 

(3.28)  R,(x)  -  O.W <»-»«’»)  . 

J  (l+a0x)(l+00'x) 


Since 


(3.29) 
we  should  prescribe 


lim  |R,(x)|  -  0/a  , 
x-*+“ 


(3.30)  a  a  0 


to  have,  from  (3.26),  (3.27),  the  stability  of  scheme  (3.9)-(3.12)  for 
the  large  eigenvalues  of  A.  We  discuss  now  the  accuracy  of  scheme  (3.9) 
(3.12)  ;  we  can  show  that 

2 

(3.31)  R3(x)  -  1-x  +  ^-{1  +  (02-a2)(262-40+l)}  +  x2  n(x), 

with  lim  n(x)  -  0.  It  follows  from  (3.31)  that  scheme  (3 . 9) — (3 . 12) 
x  -*■  0 

is  second  order  accurate  if,  either 

(3.32)  a  -  0  (-  1/2  from  (3.13))  , 
or 


«  . 
•V 
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(3.33)  e  -  1-  /2/2  -  .29289.,.,; 

scheme  (3.9)-(3.12)  is  only  first  order  accurate  if  neither  (3.32),  nor 

(3.33)  holds. 

If  one  takes  a  ■  8  ■  1/2  it  follows  from  (3.27)  that  scheme  (3.9)- 
(3.12)  is  unconditionnaly  stable  V0  e  ]0,l/2[  ;  however  since  (from 
(3.29))  we  have 

(3.34)  lim  | R. (x) |  -  1 
x  -►  +» 

the  remark  done  for  scheme  (3.3)-(3.5),  about  the  integration  of  stiff 
systems,  still  holds.  In  general,  we  shall  choose  a  and  8  in  order  to 
have  the  same  matrix  for  all  the  partial  steps  of  the  integration  me¬ 
thod,  i.e.  a, 6, 0  have  to  satisfy 

(3.35)  a0  -  8(1-20) 
which  implies 

(3.36)  a  -  (1-29)/ (1-0) ,  6  -  8/0-6). 

Combining  (3.30),  (3.36)  we  obtain 

(3.37)  0<8i  1/3. 

For  6  ■  1/3,  (3.36)  implies  a  ■  0  •  1/2. 


I 


If  0  <  0  <  1/3  and  if  a  and  8  are  given  by  (3.36)  we  have  then 


I 

i 


I  •  , 

.*  /  *.*  *.*,  «*  •*  •*,  * \  /,  •*  sm  •*,  «*,  A  •  •  •  •  •*.  A  *•  •  A  • » •  «,*  • « ’  •  » A  \A  Ai  i 


(3.38) 


lim 

n  -*  +  » 


R3(x) 


_  e  _  9 

a 


<  1 


Actually  we  can  prove  that  0  e  J  0,1/3  ]  and  a  and  8  given  by  (3.36), 
imply  the  unconditional  stability  of  scheme  (3.9)— (3. 12)  ;  moreover  if 
0  e  ]  0,1/3  [  (with  a,  8  still  given  by  (3.36)),  property  (3.38)  makes 
that  scheme  (3.9)-(3.12)  has  good  asymptotic  properties  as  n  ■*  ♦  ®  and 
for  example  is  well  suited  to  compute  steady  state  solutions. 


[f  8  -  1  -  /2/2  (resp.0  »  .25),  we  have  a«  2-/2,  8  ■  /2-1,  8/a  ■  1//2 
(reap,  a  ■  2/3,  8  ■  1/3,  8/a  ■  1/2). 

3.3.  Application  to  the  aolution  of  the  time  dependent  Navier-Stokes 


equations. 

We  discuss  now  the  application  of  the  schemes  described  in  Sec.  3.1 
to  the  solution  of  the  time  dependent  Navier-Stokes  equations  (2.1),  (2.2) 
with  the  initial  value  condition  (2.5)  ;  we  suppose  for  simplicity  that 
the  boundary  conditions  are  of  the  Dirichlet  type,  i.e. 


(3.39) 


u  ■  g  on  T  ( with 


ith  |  g.n  dr  •  0) . 


3.3.1.  A  first  operator  splitting  method. 

This  method,  which  is  directly  derived  from  the  Peaoeman-Rachford 
8oheme  (3 . 3) — (3 . 5)  is  described  as  follows  : 


(3.40) 


u  -  u  , 

•W 


-k  |k  1  1  /o  qX  I  /O 

then  for  n  a  0  and  starting  from  u  we  compute  {u  '  ,  p  '}  and 


o+J  ,  -  . 

u  by  solving 


(3.41) 


n+1/2  _  n 

~  ~  v  A  n+1/2  .  „  n+1/2  ..n+1/2  .  v  A  n  ,  n  n\  n 

— — — —  -  t  iu  +vp  -f  +  •=■  Au  -(u  .v)u 

At/2  2  ~  ~  ~  2  -  ~  ~  ~ 

„  n+1/2  n  n 
V.u  -  0  in  ft  , 


n+1/2  n+1/2  _ 

u  -  g  on  T  , 


(3.42) 


■i!.*1-  h!1 —  .  V  iu"*1.<u”*'.v>untl  -  f”*1.  i  iu"tl/2-  Vp"*1'2 

At/2  2 . 2  ~ 

n+ 1  n+ 1  _ 

u  ■  g  on  T  t 


respectively . 


3.3.2.  A  second  operator  splitting  method 


This  method  is  derived  from  scheme  (3.9)-(3.12)  and  is  described 


as  follows  : 


t.'.V 
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gr«ater  complexity  scheme  (3.43)-(3.46)  is  (per  step)  almost  as  econo¬ 
mical  to  use  as  scheme  (3. 40) -(3. 42)  ;  this  is  mainly  due  to  the  fact 
that  the  "quasi"  steady  Stokes  problems  (3.41)  and  (3.44),  (3.46) 
(actually  convenient  finite  element  approximation  of  them)  can  be  solved 
by  qui4"*  efficient  solvers  so  that  most  of  the  computer  time  used  to 
solve  a  full  step  is  in  fact  used  to  solve  the  nonlinear  subproblem. 

The  good  choice  for  a  and  6  is  given  by  (3.36)  if  one  uses  scheme  (3.43) 

(3.46)  ;  with  such  a  choice  many  computer  subprograms  can  be  used  for 
both  the  linear  and  nonlinear  subproblems,  resulting  therefore  in  quite 
substantial  core  memory  savings. 


3.4.  Application  to  eigenvalue  calculations. 

3.4.1.  Generalities.  Synopsis. 

The  main  goal  of  this  section  is  to  show  that  the  concepts  intro¬ 
duced  in  Sec.  3.1  apply  also  to  eigenvalue  calculations  at  least  for 
eyrmetric  matrices  (or  operator) .  Since  the  resulting  methods  belong  to 
the  class  of  the  so-called  inverse  'power  methods,  the  new  approach  brings 
very  little  to  the  linear  eigenvalue  problem  (for  which  a  basic  reference 
is  [23]  ),  but  it  is  nicely  suited  to  solve  nonlinear  eigenvalue  problems 
like  for  example  the  Hartree  equation  in  Quantum  Physics  (see  [24]  ,  [25] 
for  more  details) . 

3.4.2.  Formulation  of  the  problem. 

Let  A  be  a  symmetric  N  *  N,  real  matrix  ;  we  denote  by  X.  the 
eigenvalues  of A  and  we  suppose  that  X.  £  X 

We  concentrate  on  the  calculation  of  the  smallest  eigenvalue  of  A 
(i.e.  X j) .  it  is  well  known  that  Xj  satisfies 

(3.47)  X.  ■  Min  (Av,v), 

1  V£  S 

where 

(3.48)  S  ■  {v jy  e  1RN,  jjv||  -  1}  , 

N  N 

(v,w)  -  l  w^  Vv  -  w 

i*  1 

llyll  -  <v.v)1/2  . 


xwiJi-l 


(3.49) 
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Let  introduce  now  the  functionel  Ig  defined  by 

f  0  if  v  «  S, 

(3.50)  I.(v)  -  « 

S  *  /  +  »  if  v  e  S  ; 

Ig  is  the  indicator  functional  of  S.  We  clearly  have 

(3.51)  A./2  -  Min  {-(Av,v)  +  I,,(v)}  . 

1  v  e!RN  2  ~  ~  S  ~ 

Suppose  that  Ig  is  differentiable  (which  is  definitely  not  the  case) 
and  denote  by  31.  its  differential  ;  if  u  is  a  minimizer  in  problems 

d 

(3.47),  (3.51),  i.e.  an  eigenvector  of  norm  one,  associated  to  Aj,  u 
satisfies 

(3.52)  Au  +  dig (u)  -  0  . 

It  is  then  quite  natural  to  associate  to  the  nonlinear  "equation"  (3.52) 
the  initial  value  problem  below 

!du 

dT  +  +  3IS(^  "  2  * 

u(0)  -  u  , 

— 

and  to  look  for  the  eteady  state  solutions  of  (3.53),  i.e.  to  lim  u(t), 

t-*  ♦  »  ~ 

if  such  a  limit  exists.  From  the  special  form  of  (3.53),  it  is  tempting 
to  solve  it  using  the  operator  splitting  methods  discussed  in  Secs.  3.1, 
3.2.  The  resulting  algorithms  are  described  in  the  following  Sec.  3.4.3. 

3.4.3.  Solution  of  (3.52),  via  (3.53),  by  operator  splitting  methods. 
We  apply  now  the  operator  splitting  methods  of  Sec.  3.1  with 
Aj  ■  3Ig,  A£  •  A  and  f  ■  0. 

Application  of  the  Peaceman-Rachford  scheme  (3.3)-(3.5)  : 

(3.54)  u°  “  uQ, 


then  for  n  *  0  with  un  known  we  compute  successively  un+ly/2  and  un+1  by 
n+l/2_  n 

<3-55)  *  5is<h”*1/2)  *  *  5°  ■  $  . 

n+l_  n+1/2 

(3.56)  At/2 -  +  Slg(un+1/2)  ♦  A  un+1  -  0  . 

Application  of  the  Douglas-Rachf ord  scheme  (3.6)-(3.8)  s 

(3.57)  u°  -  u  , 

then  for  n  a  0,  with  un  known  we  compute  successively  Qn+*  and  un+1  by 
Gn+1-  un 

(3.58)  ~  Jt'~-  +  3lg(fin+1)  +  A  un  -  0, 

n+l_  n 

(3.59)  ♦  3I_ (fln+1)  +  A  un+1  ■  0  . 

ilk  0  w  m*  m* 

Application  of  the  three  stages  scheme  (3.9)-(3.12)  : 

(3.60)  u°  -  u  , 

-O 

then  for  n  i.  0  ,  with  un  known  we  compute  successively  un+®,  un+1  ® 
and  un+ 1  by 

n+9_  n 

(3.61)  -  6At'  -  +  3ls(un+9)  +  A  un  -  0, 

n+1-0  _  n+0 

(3.62)  ^ - = -  +  31. (un+0)  +  A  un+1"6  -  0, 

(1-26)  At  b  ~ 

n+1  n+1-0 

U  —  u 

(3.63)  ■  —  . —  +  3l_(un+S  +  A  un+1  ■  0  . 

6At  S  - 

In  order  to  derive  more  practical  formulations  of  the  above  algorithms, 


let  consider  algorithm  (3.54)-(3.56)  (the  conclusion  to  be  obtained  will 
hold  also  for  the  two  other  methods)  : 


Y , ■  ,%  ,%  ,*•  . %  ,n V . ,  Z'. \* ,  •.*  •  .  . .  t  v  *>  v  v  v  v  .*,*.* 


KMMtUa 


We  observe  that  (3.55)  is  in  fact  a  necessary  condition  of  optima¬ 
lity,  for  the  following  minimization  problem 


(3.64) 


Find  un+1/2  e  s  euah  that 


Jn (°n+1/2)  *  Jn<v),  Vv  €  S, 
u  **  n 


Jn<v)  -  ±||v||  2  -(u11-  A  u“,  v)  . 

Since  |jv||  •  1,  Vv  e  S,  the  solution  of  (3.64)  is  given  by 

un+l/2  ~Q“  6  li*1 

-  'll  u°-  4  A  u”||  • 

On  the  other  hand  it  follows  from  (3.55)  that 

At  ,  n+1/2.  n  At  .  n  n+1/2 

—  ®Ig(“  )  *  u - 2  £  ~  -  u  . 

which,  combined  to  (3.56),  implies 


(I  *  4  A)u"*‘  .  2u"*1/2  -  („“  -  ^  *  „»). 

<V  4  **»  /  At 


Collecting  the  above  results  we  obtain  the  following  practical  formulation 
for  algorithm  (3.54)-(3.56) . 

Practical  formulation  of  algorithm  (3.54)-(3.56)  : 


(3.65) 


u  ■  u  , 
~  ~o 


then  for  n  a  0,  compute  pn  and  un+* ,  from  un,  by 


(3.66) 


(3.67) 


n  n  At  .  n 
p  •  u - 7  A  u  , 

^  -W  4  -  A, 


»”'  -  a  ♦  a>-1 c — 2__  -  i)  p°  . 

2  '  II  pn  II 
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Usiog  similar  calculations  we  should  obtain  the  following  formulations 
for  (3.57)-(3.59)  and  (3.60)-(3.63)  : 

Practical  formulation  of  algorithm  (3.57)-(3.59)  s 

(3.68)  u°  -  u  , 

-  -O 

then  for  n  z  0,  compute  pn  and  un+^,  from  un,  by 

(3.69)  pn  -  un  -  At  A  un  , 

(3.70)  uwl  -  (I  ♦  At  A)"1  {un  ♦  (-i-  -  l)p°)  . 

*  Up  ll  ' 

Practical  formulation  of  algorithm  (3.60)-(3.63)  s 

(3.71)  u°  -  u  . 

•*»  —  o 

^  n  n+0  n+1-6  n+1-0  n+1  »  n  , 

then  for  n  i  0,  compute  p  ,  u  ,  u  ,  p  ,  u  ,  from  u  by 


(3.72) 

(3.73) 

(3.74) 


p°  *  un  -  0At  A  un  , 


n+6 


u 


l!pnll 


n+1-0  /T  ^  1 1  OONA-  A N ■ 1 / 1-6  n+0 

u  -  (I  +  (1-20)  At  A)  (— g—  u 

-  (I  ♦  d-20)At  A) ~ 1  (—jp  — — 

Ip  I 


1-20  n, 

P  )  - 


1-20,  n 
— >?  * 


(3.75) 

(3.76) 


n+1-0  n+1-0  nA.  .  n+1-0 

p  ■  u  -  0At  A  u  , 


n+1 


n+1-0 

P 

n+1-0 

P  I 


From  the  above  relations  the  three  algorithms  which  have  been  described 
appear  as  variations  of  the  well-known  inverse  power  methods  (see,  e.g., 
[233  for  a  detailed  analysis  of  such  methods).  From  our  numerical  expe¬ 
riments  the  various  vector  sequences  generated  by  the  above  algorithms 
converge  linearly  to  an  eigenvector  associated  to  the  smallest  eigenvalue 
Xj,  while  (A  un,  un)  converges  quadratically  to  Xj.  Actually  the  conver- 


« v  *  „  .  ...  1 

r*  .*■ ,  . •  v *,•  */ v  * '«* *.*  v *•  •*. «*. ■*.  •*,<*. •*. ^ . •” , •*« 


gence  is  still  good  if  Aj  is  a  multiple  eigenvalue.  It  is  interesting 
to  observe  that  the  fastest  algorithm  is  (3.73)— (3. 76) ,  then  (3.65)- 
(3.67),  then  finally  (3.68)— (3.70)  ;  these  results  agree  with  the 
analysis  done  for  the  trivial  linear  problem  that  we  considered  in 
Sec.  3.2. 

In  the  following  sections  we  go  back  to  the  Navier-Stokes  equations 
and  their  numerical  treatment. 


4.  Least  squares  coniusate  gradient  solution  of  the  nonlinear  subproblems 


obtained  from  the  time  discretization  of  the  Navier-Stokes  E 


the  methods  of  Section  3.3. 


4.1.  Classical  and  variational  formulation.  Synopsis 


At  each  full  step  of  the  operator  splitting  methods  (3.40)-(3.42) 
and  (3.43)-(3.46)  we  have  to  solve  a  nonlinear  elliptic  system  of  the 
following  type 


(4.1) 


au  -  vAu  ♦  (u.V)u  ■  f  in  Si, 


u  -  g  on  T, 


where  a  and  v  are  two  positive  constants  and  where  f  and  g  are  two 
given  functions  defined  on  SI  and  T  ,  respectively.  We  do  not  discuss 
here  the  existence  and  uniqueness  of  solutions  for  problem  (4.1).' 

We  introduce  now  the  following  functional  spaces  of  Sobolev’s  type 


(4.2)  H1  (ft)  -  £  L202),  e  L2  («) ,  Vi 

l 

(4.3)  Hq(S1)  -  {<j>|<{>  e  Hl(fl),  $  -  0  on  H  , 


1.....N}  , 


(4.4)  Vq  -  (H^(fl))N  , 


(4.5) 


V  -  {v|v  e  (HJ(fi))N,  v  -  g  on  T)  ; 


if  g  is  sufficiently  smooth  then  V  is  nonempty. 

**  S 

We  shall  use  the  following  notation 


and  if  u  -  v  -  ,  then 


u.v  ■  Z  u.v., 

*v  <v  •  «  11 

1*1 

N  N  N  3u.  &v. 

Vu.Vv  —  Z  Vu. .  Vv.  “  Z  Z  T--  'g"«  . 

*•"*  .  •  **  1  1  i«i.  0X«  vX« 

i-i  i-i  j-i  j  j 

Using  Green's  formula  we  can  prove  that  for  sufficiently  smooth  functions 

u  and  v  belonging  to  (H*(ft))N  and  V  ,  respectively,  we  have 
~  ~  o 


►  - 6)  -  [  Au.v  dx  -  I  Vu.Vv  dx  . 


It  can  also  be  proved  that  if  u  is  a  solution  of  (4.1),  belonging  to  V  , 

O 

it  is  also  a  solution  of  the  following  nonlinear  variational  ■problem 


(4.7) 


u,  Vg, 


a  f  u.v  dx  +  v  [  Vu.Vv  dx  +  (  ((u.V)u).v  dx  -  [  f.v  dx  ,VvcV 
J  ~  J  ~~  ~~  J  ~  ~  ~  ~  ;  -  ~  ~  c 


and  conversely.  We  observe  that  (4.1), (4. 7)  is  not  equivalent  to  a 
problem  of  the  Calculus  of  Variations  since  there  is  no  functional  of 
v  with  (v.V)v  as  differential  ;  however  using  a  convenient  least  squares 
formulation  we  shall  be  able  to  solve  (4.1),  (4.7)  by  iterative  methods 
originating  from  Nonlinear  Programming,  such  as  conjugate  gradient  for 
example. 

4.2.  Least  squares  formulation  of  (4.1),  (4.7). 


Let  v  c  V  ;  from  v  we  define  y  (-  y(v))  e  V  as  the  solution  of 
~  g  ~  ~  ~  ~  o 


(4.8) 


ay  -  vAy  -  av  -  vAv  +  (v.V)v  -  f  in  ft  , 


y  -  0  on  T. 


We  observe  that  y  is  obtained  from  v  via  the  solution  of  N  uncoupled 
linear  Poisson  problems  (one  for  each  component  of  y),  using  (4.6)  it 
can  be  shown  that  problem  (4.8)  is  equivalent  to  the  linear  variational 
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f  Find  y  e  Vq  such  that ,  Vz  e  Vq,  we  have 
(4.9)  I  a  I  y.z  dx  +  v  I  Vy.Vz  dx  ■  a  I  v.z  dx  +  v  I  Vv.Vz  dx 

J  J  ~  J  «w<w  J  •»  «-  I  *s»  «w 


)  n  n  a  n 

|  +|  ((v.V)v).z  dx  -  j  f.z  dx, 

V.  «  n 

which  has  a  unique  solution.  Suppose  now  chat  v  is  a  solution  of  the 
nonlinear  problem  (4.1),  (4.7)  ;  the  corresponding  y  (obtained  from 
the  solution  of  (4.8),  (4.9))  is  clearly  y  •  0  ;  from  this  observation 
it  is  quite  natural  to  introduce  the  following  ( nonlinear )  least  squares 
formulation  of  (4.1),  (4.7) 


(4.10) 


Find  u  e  V  such  that 
-  8 


J(u)  S  J(v),  ¥v  fi  V  . 


where  J  :  (H1  (12)  )N  -*■  ]R  is  the  function  of  v  defined  by 


(4.11) 


J(v)  -  -r  I  (alyl2  ♦  vJVy)2}  dx 


where  y  is  defined  from  v  by  (4.8),  (4.9).  We  observe  that  if  u  is  solu- 
tion  of  (4.1),  (4.7),  then  it  is  also  a  solution  of  (4.10)  such  that 
J (u)  ■  0  j  conversely,  if  u  is  a  solution  of  (4.10)  such  that  J(u)  ■  0, 
then  it  is  also  a  solution  of  (4.1),  (4.7). 

4.3.  Conjugate  gradient  solution  of  the  least  squares  problem  (4.10). 
4.3.1.  Description  of  the  algorithm. 

We  use  the  Pola’k-Ribidre  version  (see  [26]  )  of  the  conjugate 
gradient  method  to  solve  the  minimization  problem  (4.10)  ;  we  have  then 
(with  J’ (v)  the  differential  of  J  at  v) 

Step  0  :  Initialization 


(4.12) 


u  e  V  ,  given  ; 


we  define  then  g°  ,  w°  £  VQ  by 


(4.13) 


§°  6  V 


a  [  g°.z  dx  +  v  f  Vg°,Va  dx  ■  <J'(u°),z>  ,  Vz  e  V  , 

j  j  %.  •»  I  ««  «W  «V  O 

(,  a  n 


(4.14) 


o  o 

Jf  "  I  » 


respectively .  □ 

Tfcen  /or  n  *  0,  assuming  that  un,  wn,  gn  are  known,  we  obtain  un+1, 

n+1  n+1  , 

g  .  w  by 


Step  1  :  Descant 


(4.15) 


(4.16) 


Find  X  e  K  » 
n 


J(un-  X  wn)  a  J (un-X  w11),  VX  (  K  , 

•«  Q  %  —  -w 


n+1  n  .  n 
u  -  u  -  X  w 

—  Q 


Step  2  J  Calculation  of  the  new  descent  direction 


(4.17) 


(4.18) 


Find  gn+1  «  V  .  such  that 
«*»  o 


dx  »  <J,(un+S»z>  *  Vz  c  V  , 

-*»  ••  •»  o 


a  [  gn+1.z  dx  +  v  I  Vgn+1.Vz 

n  n 

a  [  g°+1 • (gn+1-gn)dx  +  v  [  Vgn+l.V(gn+1-gn)dx 


Yn" 


a  [  |gn| 2  dx  +  v  [  |Vgn|2 


dx 


(4.19) 


n+1  n+1  n 

w  •  g  +  y  w 
—  —  n  - 


□ 


Do  n  -  n+1,  go  to  (4.15). 


As  we  shall  see,  applying  algorithm  (4. 12)— (4 . 19) ,  to  solve  the 
least  squares  problem  (4.10),  requires  the  solution  at  each  iteration 
o£  exactly  three  Diriahlet  systems  (i.e.  3N  scalar  Dirichlet  problems) 
associated  to  the  elliptic  operator  al-vA  . 


V  V  V  V.V,  •V'"- *  .<  •  ’  •  ■  o  . > .  •- 


4.3.2.  Calculation  of  J*. 

A  most  important  step  when  making  use  of  algorithm  (4. 1 2>  —  (4 . 19) 
to  solve  problem  (4.10),  is  the  calculation  of  <J'(un+1),z>  at  each 
iteration  ;  we  should  easily  prove  (see,  e.g.,  [5]  ,  [22])  that  J' (v) 
can  be  identified  with  the  linear  functional  from  Vq  to  H,  defined  by 

(  <J'(v),z>  ■  a  y.z  dx  +  v  f  Vy.Vz  dx 

(4.21)  )  ~  ~  J°“"  j*~~~ 

I  +  |  y. (v.V)z  dx  +  |  y.(z.V)v  dx,  ¥z  e 

where  0  ft 

N  3w. 

(v.V)w  -{  Z  v.  T-i  ,  ,  Vv,  w  ; 

j-1  J  dxj  'i-1  ~  ~ 


<J,<Y),5>  ^as  therefore  a  purely  integral  representation  which  is 
of  major  importance  in  view  of  finite  element  (or  spectral )  implemen¬ 
tations  of  algorithm  (4. 12)-(4. 19) .  From  the  above  results,  to  obtain 
<J'(un+1),z>  we  should  proceed  as  follows  : 

(i)  Compute  yn+1,  associated  to  un+I  by  (4.8),  (4.9),  as  indicated 

*** 

in  Sec.  4.3.3,  below. 

(ii)  We  obtain  then  <J' (un  ^),z>  by  taking  v  ■  un+^  and  y  m  yn+^ 
in  (4.21). 

4.3.3.  Calculation  of  Xq.  Further  comments  on  algorithm  (4. 12) —(4 . 19) . 
A  problem  of  practical  importance  is  the  calculation  of  X  .  Let 

XI 

denote  by  y  (X)  the  solution  of  (4.8),  (4.9)  associated  to  v  -  un-  Xw11  ; 
we  clearly  have 


(4.22)  ya(0)  -  y11,  yn(X  )  -  yn+1  , 

*V  •»  U 


and  also 


(4.23)  yn(X)  -  yn  -  Xya  ♦  X2  ya  , 


where  ya,  ya  are  the  solutions  of 


V  *  V*” 

v.  .  -v'v;. . . . ; 

»  .  •  a*  a »  -  -  k.  *  -  •  .  •  k  •  fc  *  k.  f .  LT*  .Vl  >  k  ^  m  • 

:•  .-.v; 

»*>  \*»V»  , 

aU  v  v'  a’  o 

S 

<•«  w  n  n 

«*rvr.T 

■V 

-23- 


(4.24) 


ay”-  vAy”  -  aw”  -  vAw”+  (u”.V)w”+(w”.V)un  in  Q, 

■V  1  —1  -«»  «*  •*  **•  *  ■»  —  •*» 


y}  »  0  on  T, 


(4.25) 


ay”  -  vAy”  ■  (w”.V)w”  in  !2 , 
— z  —  — 


y2  -  0  on  T, 


respectively.  Since 


(4.26) 


J(u”-  Xw”)  -  \  J  {ajy”(X)|2  +  v|Vy”(X)|2}  dx  , 


n...  1 2 


the  function  X  J(u”-  Xw”)  is,  from  (4.23X  a  quartio  polynomial  in  X 

•«» 

that  we  shall  denote  by  j  (X)  ;  X  is  therefore  solution  of  the  cubic 

n  n 

equation 


(4.27)  j^(X)  -  0  . 


He  shall  use  the  standard  Newton  method  to  compute  Xq  from  (4.27), 
starting  from  X  -  0.  The  resulting  algorithm  is  given  by 


(4.28) 


X°  -  0, 


lr+1  Ir 

then  for  k  2  0,  we  obtavn  X  from  X  by 


(4.29) 


, k+1  . k 

X  ■  X  - 


ij(xk)  ' 


In  our  calculations,  we  always  observed  a  very  fast  convergence  of 
algorithm  (4.28),  (4.29).  Once  X  is  known,  we  know  y”+*  since  (from 
(4.22))  y”+1  -  y”(X  ) . 

If  we  count  now  the  number  of  Dirichlet  systems  for  al-vA,  to  be  solved 
at  each  iteration,  we  observe  that  we  have  to  solve  only  three  such 
systems,  namely  (4.24),  (4.25)  and  then  (4.17)  (to  obtain  g”+*)  ;  this 
number  is  optimal  for  a  nonlinear  problem  since  the  solution  of  a  linear 


problem,  by  least  squares-preconditionad  conjugate  gradient,  requires 
the  solution  at  each  iteration  of  two  linear  systems  associated  to 
the  preconditioning  operator. 

From  the  above  remarks,  it  appears  clearly  that  the  practical 
implementation  of  algorithm  (4. 1 2)  —  (4 . 19)  will  requires  an  efficient 
(direct  or  iterative)  elliptic  solver,  like  one  of  those  discussed 
in  [27]  ,  [28]  .  As  a  final  comment,  we  would  like  to  mention  that 
algorithm  (4. 12)-(4. 19)  (in  fact  its  finite  dimensional  variants)  is 
quite  efficient  ;  when  used  in  combination  with  the  operator  splitting 
methods  of  Sec.  3,  three  to  five  iterations  suffice  to  reduce  the  value 
of  the  cost  function  by  a  factor  of  10^  to  10^  ;  however  in  view  of 
other  applications  we  are  testing  now  some  of  those  methods  combining 
the  features  of  conjugate  gradients  and  quasi-Newton  algorithms,  such 
as  the  methods  discussed  in  [29],  [300  (the  results  recently  obtained 
for  the  calculation  of  transonic  potential  flows  containing  shocks ,  look 
very  promissing) . 

5.  Solution  of  the  "quasi"  Stokes  linear  subproblems. 

5.1.  Generalities.  Synopsis. 

At  each  full  step  of  the  splitting  methods  discussed  in  Sec.  3.3 
we  have  to  solve  one  or  two  linear  problems  of  the  following  type 

r  otu  -  vAu  +  Vp  ■  f  in  ft, 

mm  mm  mm  mm 

(5.1)  V.u  ■  0  in  ft, 

^  u  -  g  on  T  ( with  |  g.n  dr  -  0), 

where  a  and  v  are  two  positive  constants,  and  where  f  and  g  are  two 

given  functions  defined  on  (1  and  T  ,  respectively.  We  recall  that  if 

f  and  g  are  sufficiently  smooth,  then  problem  (5.1)  has  a  unique 

solution  in  V  x  (L2(fl)/K)  (with  V  still  defined  by  (4.5)  ;  p  e  L2(«)/ * 

o  g 

means  that  p  is  defined  only  to  within  an  arbitrary  constant).  We  shall 
describe  below  two  iterative  methods  for  solving  (5.1),  quite  easy  to 
implement  using  finite  element  or  spectral  methods  (more  details  are  given 


[5]  ,  [22]  ,  together  with  convergence  proof*  ;  more  methods  ere 
discussed  in  [5  3  ) . 


5.2.  A  conjugate  gradient  algorithm  for  solving  (5.1). 

A  complete  justification  of  the  following  algorithm  is  given  in 
[5  ]  ;  we  can  say  briefly  that  eliminating  u  in  (5.1),  it  appears  that 

•m 

p  is  solution  of  a  linear  functional  equation  associated  to  an  opera¬ 
tor  which  is  self  adjoint  and  strongly  elliptic  from  B  onto  H,  where 

H  -{q|q  e  L2(fi),  j  q  dz  ■  0}  . 

n 

Such  properties  justify  a  conjugate  gradient  solution  of  (5.1)  and 
lead  to  the  algorithm  below 


(5.2)  p°  e  L2(fl),  given. 


(5.3) 


au°  -  vAu°  »  f  -  7p°, 

•*» 

u°  ■  g  on  T, 


(5.  A) 


O  _  _  o 
8  •  V.u  , 


(5.5) 


O  0 

w  ■  g 


-  _nnn..,  ^n+ln+ln+1 

Then  for  n  *  0,  p  ,  g  ,  w  bexng  hxown,  compute  p  ,  g  ,  w  ae 

follows  : 

Solve 


(5.6) 


axn  -  vAx°  -  -Vwn, 
Xn  -  0  on  r, 


(5.7) 


I.' 


gn|2  d* 


P_  ■ 


n  fn.ni 
V.x  w 
- 


dx 


n+ 1  n  _  n 

p  -  p  -  pQ  w 


(5.8) 


» 


V  V 


sniBSB 


(5‘5)  Y"  ‘  I  uni2  <*  ’ 


(5.10)  wn+l  -  gn+1  +  Yn  ^  • 


Do  n  ■  n+1,  to  (5.6). 

Ue  observe  that  algorithm  (5.2)-(5. 10)  requires  the  solution  at  each 
iteration  of  only  one  Dirichlet  system,  namely  (5.6).  It  can  be  proved 


(5.11)  lim  {u  ,  pn}  -  {u,  p  } 

n  -+•  +00 


where  {u,  p  }  is  that  solution  of  (5.1)  such  that 
•»  o 

J  PQ  dx  "  J  P°  d*  • 

ft  ft 

5.3.  A  second  iterative  method  for  solving  (5.1). 

This  method  is  defined  as  follows  (with  r  a  nonnegative  parameter)  : 

(5.12)  p°  £  L2(ft),  given , 

then  for  n  i  0  ,  define  un+1  and  pn+1  ,  from  pn  by 


(5.13) 


(5.14) 


au11  -  vAun  -  r  V(V.un)  ■  f  -  Vpn, 


n  r 

u  ■  g  on  T 


n+1  n  „  „  n 
p  -  p  -  p  v.u  . 


The  above  method  is  related  to  the  artificial  compressibility  methods  of 
Chorin  and  Yanenko  since  (5.14)  can  be  considered  as  a  discretization  of 


♦  V.u  -  ° 


(  p  then  plays  Che  role  of  a  cine  seep) . 

About  Che  convergence  of  (5. 12)-(5. 14)  we  should  prove  (see,  e.g., 
[5,  Chapcer  7  ]  for  such  a  proof)  Che  following 


Proposicion  5.1.  :  Suppose  that 

(5.15)  0  <  p  <  2(r+  £)  ; 
we  have  then 

(5.16)  lim  {un,  pn)  ■  {u,  p  }  strongly  in  (H* (ft) )^  x  L2(ft) 

*V  •  •*  O 

n  -*  +  » 

where  {u,  p  }  is  the  solution  of  (5.1)  such  that 

•V  O 

}  po  dx  “  J  P°  • 

ft  ft 

Moreover  the  convergence  is  linear  (i.e.  ||un-  u||  cmd 

||pn-  p  ||  ,  converges  to  zero  as  fast,  at  least,  as  geometric 
IT  (ft) 

sequences ) . 

Remark  5.1: ( About  the  choice  of  p  and  r  )  :  We  should  use  p  ■  r  in 

pracCice,  since  ic  can  be  proved  in  ChaC  case  ChaC  Che  convergence 

racio  of  algorichm  <5.12)— (5. 14)  is  0(r  *),  for  large  values  of  r.  In 

2  4 

mosC  applicacions  caking  r  -  10  v  Co  10  v  we  have  a  practical  conver¬ 
gence  of  algorichm  (5. 12) — (5 . 14)  in  3  Co  4  iCeracions.  There  is  however 
a  praccical  upper  bound  for  r,  since,  for  Coo  large  values  of  r,  problem 
(5.13)  will  be  ill-conditioned  and  ics  praccical  solucion  sensicive  co 
round-off  errors. 

Remark  5.2  :  If  r  -  0,  problem  (5.13)  reduces  Co  Che  solucion  of  N  un¬ 
coupled  (one  for  each  componenc  of  un)  scalar  Dirichlet  problems  for 
al-vA  ;  if  r  >  0  Che  N  components  of  un  are  coupled  by  V(V.un),  making 

•w  A.  % 

Che  solucion  of  (5.13)  much  more  cosCly.  In  face,  Che  ellipCic  operacor 


in  the  left  hand  side  of  (5.13)  is  very  close  to  the  linear  elasticity 
operator,  and  close  variants  of  it  occur  naturally  in  compressible  and/or 
turbulent  viscous  flow  problems. 


Remark  5.3  :  Other  methods  for  solving  (5.1)  are  discussed  in  [5  ]  ,  [31  3, 
[32  1 


6.  Finite  Element  Approximation  of  the  Time  Dependent  Navier-Stokes 
Equations. 

We  shall  describe  in  this  section  a  specific  class  of  finite 
element  approximations  for  the  time  dependent  Navier-Stokes. 

Actually  these  methods  which  lead  to  continuous  approximations  for 
both  pressure  and  velocity  are  fairly  simple,  and  some  of  them  have 
been  known  for  years.  They  have  been  advocated  for  example  by  Hood 
and  Taylor  (see  [33]  ).  Other  finite  element  approximations  of  the  in¬ 
compressible  Navier-Stokes  equations  can  be  found  in  [1],  [2],  [41  [5], 
[34  ]  (see  also  the  references  therein) . 


6.1.  Basie  hypotheses.  Fundamental  discrete  spaces. 

2 

We  suppose  that  (1  is  a  bounded  polygonal  domain  of  K  .  With 
•ig’k  a  standard  finite  element  triangulation  of  £1  ,  and  h  the  maximal 

length  of  the  edges  of  the  triangles  of  we  introduce  the  following 

discrete  spaces  (with  P^  ■  space  of  the  polynomials  in  two  variables  of 
degree  £  k)  : 


(6.!)  h;  -  {qh|qh  €  C° (S2)  ,  qJT  €  P,,  VT  tfh  }  , 

(6.2)  Vh  -  {vh|vh  e  C°(5)  x  C°(i5),  yh|T  t  Pj  x  Pj,  ?I  cfh  ), 

(6*3)  Voh  “{?h  6  V  *h  "  2  on  F}  "  Vh  n  Vo  * 

Two  useful  variants  of  V,  (and  V  ,  )  are  obtained  as  follows  s  either 

h  oh 


(6.4)  Vh  -  {yh|yh  c  C°(fl)  x  C°(£2),  vJT  £  Pj  xPj,  VTt^}, 


or  (this  space  has  bean  introduced  in  [35  ]  ) 


<6-5>  vh  •  C°<S)  *  c°<s>'  vhlT  «  P?T  *  p*r  w  '^h}  • 

In  (6.4),  c(fh  is  this  triangulation  of  ft  obtained  from  9?^  by  joining 
the  midpoints  of  the  edges  of  I  c  ^  ,  as  shown  on  Fig.  6.)  ;  we  have 
the  same  global  number  of  unknowns  if  we  use  defined  by  either  (6.2) 
or  (6.4),  however  the  matrices  encountered  in  the  second  case  are  more 
compact  and  sparse.  In  (6.5),  P*T  is  the  subspace  of  P^  defined  as  fol¬ 
lows 

(6.6)  P*T  -  {q|q  -  qj  ♦  X$T,  with  qjC  Pj.XeK,  and  <|»T  €  Pj, 

♦j  ■  0  on  3T,  4>t(Gt)  -  1} 

where,  in  (6.6)  is  the  oenteotd  of  T  (see  Fig.  6.2  below).  A  function 
like  4>t  is  usually  called  a  bubble- function. 


6.2.  Approximation  of  the  boundary  conditions. 

If  the  boundary  conditions  are  defined  by 

(6.7)  u  -  g  on  T  with  g.n  dr  ■  0, 

r 

it  is  of  fundamental  importance  to  approximate  g  by  g,  such  that 
I  g^.n  df  ■  0.  The  construction  of  such  g^  is  discussed  in  [5,  Appendix  3  ] 
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6.3.  Space  approximation  of  the  time  dependent  Navier-Stokes  equations. 
Using  the  spaces  H^,  and  V  ^  we  approximate  the  time  dependent 
Navier-Stokes  equations  as  follows  : 


! 


(6.8) 


(6.9) 

(6.10) 
(6.11) 


Find  (uh(t),  Ph(t)>  e  Vh  x  h^,  Vt  k  0,  such  that 

J  ir  •  ih  *** +  v  J  SvSk**  +  j  (vPv?h dx 

n  n  n 

♦ 1  ’VS h  1111  ■  I  !h-ih dx>  ’ih  6  voh  • 


J  !-«h  \  d*  ■  °>  v»h  «  %  • 


a 


Sh  ‘  *h  °»  r  • 

?h(‘-0)  *  Soh  Mth  Soh  6  V  1 


in  (6.8)-(6. 1 1) ,  f.  ,  u  ,  and  g.  are  convenient  approximations  of  f,u  and  g, 

-on  — n  —  — o  — 

respectively. 

6.4.  Time  discretization  of  (6.8)-(6.11)  by  operator  splitting  methods. 
We  consider  now  a  fully  discrete  version  of  scheme  (3. 40) -(3. 42) 
discussed  in  Sec.  3.3  ;  it  is  defined  as  follows  (with  At  as  in  Sec.  3.3) 


(6.12) 


Sh  *  Soh  • 


then  for  >i0,  compute  {from  u“)  {u“'fl/2,  p“*1/2)«  Vh  x  H l,  and  then 
U°+1  e  vh,  by  solving 


/  n+1/2  n 

'  '  %  "  Hh 


1 

n 


At/2 


(6. 13) j  < 


•Sh dx  *  f  f  sr/2.ah d*  *  f  LC‘/2-Sh  d*  • 

n  n 


"  I  ^h+1/2^h  dx  “  i  \  !“h-^hdx  ‘  |  (&Z>Hh-lhdx’  €  voh « 
n  n  n 


(6.1 3) 2  j  V.u“+1/2  qh  dx  -  0,  Vqh  e  , 

V  ~ 


.v.V.v  ^  V.*/  «.'  N*  *.  \ 


zr  1  11+1/2  ••  _UT*  I  *  4  UTi/f  kk-r  4  /  4m _  ft 

(6.13)j  uh  £  Vh,  ph  £  V  uh  -  lh  on  r  » 


n+1/2  ,  „1  n+1/2  _  n+1/2 

ph  e  V  "h  lh 


(6.14) 


(6.15)  u 
respectively. 


n+l  n+1/2 

[  ! 

lh  ”  ~h 

— — -  .v.  dx 

+  t! 

n  n+l  , 
Vu,  • 

n 

At/2  ~h 

2  J 

— n 

n 

*  j 

z  n+l  n+l  . 

(Sh  -PSh  -ih  +1  ‘ 

.n+l 

f.  .V. 

n 

£1 

V 

“  2 

I  *;r‘/2-s!h  *■  -  j 

!C‘/2 

£1 

£1 

n+1 

„  n+l  n+l 

2h 

e  V,  ,  u,  -  g. 

n  ~n  *v  n 

on  r, 

1 

~h  ’  ~h  c  oh 


The  same  techniques  apply  to  the  space  discretization  of  scheme 
(3.43)-(3.46)  (see  [5  ]  for  more  details). 

The  solution  of  the  various  subproblems  encountered  at  each  step  of 
(6. 12)— (6 . 15)  can  be  done  by  the  discrete  variants  of  the  method  dis¬ 
cussed  in  Secs.  4  and  5,  see, again,  [5]  for  more  details. 


7.  Numerical  experiments. 

Ue  illustrate  the  numerical  methods  described  in  the  above  sections 
by  the  presentation  of  the  results  of  numerical  experiments  where  these 
methods  have  been  applied  to  simulate  some  incompressible  viscous  flows 
of  practical  interest.  All  the  calculations  which  follow  have  been  done 
using  the  finite  element  method  associated  to  the  discrete  spaces  defined 
by  (6.1),  (6.3),  (6.4)  (some  preliminary  and  promissing  results  have  been 
obtained  recently,  using  defined  bv  (6.5),  (6.6)). 

7.1.  A  first  class  of  test  problems 

We  consider  the  solution  of  the  Navier-Stokes  equations  for  the 
flow  of  an  incompressible  viscous  fluid  in  the  channel  with  a  step  of 
Figure  7.1.  We  have  selected  this  problem  since  it  is  a  quite  classical 
and  significant  test  problem  for  Navier-Stokes  solvers  (see  [34  ]  for 


a  comparison  of  various  methods  for  solving  this  test  problem) .  The 
finite  element  triangulations  used  for  these  calculations  are  shown 
on  Figure  7.1.  ;  the  coarse  (resp.  fine )  one  is  used  for  approximating 
the  pressure  (resp.  the  velocity) . 


r„  ?. 

Nodes  619  2346 

Triangles  1109  4436 

Cholesky's  coefficients  21654  154971 


Figure  7.1. 


We  have  also  indicated  on  Fig.  7.1.,  the  number  of  nodes,  triangles, 
and  nonzero  Cholesky  coefficients  of  the  matrix  approximating  -  A  (resp. 
al-vA)  on  cl^h  (resp.  .  The  methods  described  in  the  above  sections 
have  been  applied  to  compute  the  steady  state  solutions  of  (2.1),  (2.2) 
for  the  following  boundary  conditions 


f  u  satisfies  Poiseuille  velocity  profiles  at  the  entrance 
|  and  exist  of  the  channel  and  is  equal  to  0  elsewhere  on  T. 


We  have  taken  Re  ■  100  and  191.  The  numerical  results  agree  quite  well 
with  those  in  [34  ],  [36  ]  and  show  a  clear  superiority  of  the  schemes 
derived  from  (3.43)— (3.46) ,  over  the  schemes  derived  from  (3.40)-(3.42) . 
For  the  scheme  (3. 43) -(3. 46)  we  have  tested  6  -  .25  and  0  ■  1-/2/2  (for 
the  same  At)  ;  the  convergence  to  the  steady  state  is  faster  with  the 
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second  value  of  6  (for  this  class  of  problems  at  least).  We  have  shown 
on  Fig.  7.2  (resp.  Fig.  7.3)  the  stream  lines  (reap,  the  ieobar  lines) 
of  the  steady  state  solutions  corresponding  to  Re  ■  100  and  Re  •  191. 

We  observe  that  the  sine  of  the  recirculation  region  increases  with  Re. 


Stream  lines 


7.2.  A  second  class  of  test  problems 


The  second  test  problem  that  we  considered  is  much  more  complica¬ 
ted  than  the  first  one,  since  it  concerns  the  simulation  of  an  incom¬ 
pressible  viscous  flow  around  and  inside  a  (two-dimensional)  nozzle  at 
high  incidence  (40  degrees) ,  and  at  Re  ■  750  (the  characteristic  lengtn 
being  the  distance  between  the  walls  of  the  nozzle) .  We  used  the  same 
kind  of  finite  element  approximation  than  in  Sec.  7.1.  Figures  7.4,  7.5 
show  the  details  of  the  triangulations  *6*^  and  cg>^,  respectively,  close 
to  the  air  intake.  Figures  7.6-7.10  show  the  stream  lines  and  the  vortex 
patera  of  the  flow  at  t  ■  .0,  .2,  .4,  .6,  .8,  the  initial  velocity  being 
associated  to  the  corresponding  steady  Stokes  flow,  and  a  suction  phe¬ 
nomenon  being  simulated  inside  the  nozzle. 


Figure  7.4. 
Pressure  grid 
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8.  Conclusion 

We  have  presented  in  this  paper  numerical  methods  for  solving  the 
time  dependent  incompressible  Navier-Stokes  equations.  One  of  the  key 
ingredient  to  make  the  solution  process  faster ,  and  lees  demanding  from 
the  computer  storage  point  of  view,  was  definitely  the  use  of  operator 
splitting  methods  in  order  to  decouple  the  two  main  difficulties  of  the 
problem,  namely  the  nonlinearity  and  the  incompressibility.  Others  split* 
ting  methods  for  the  incompressible  Navier-Stokes  equations  are  discussed 
in  [37]  -  T40  ]  ,[  1  3,  C 1 7  ]  .  Another  important  reference  for  operator 
splitting  methods  is  [41 ]  . 

Application  of  these  splitting  methods  to  nonlinear  eigenvalue 
problems  will  be  discussed  in  a  forthcoming  paper. 
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